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ABSTRACT 

A review is presented covering both finite difference and finite element 
analysis of small amplitude (linear) sound propagation in straight and variable 
area ducts with flow, as might be found in a typical turbojet engine duct, muf- 
fler, or industrial ventilation system* Both "steady" state and transient 
theories are discussed. Emphasis is placed on the advantages and limitations 
associated with the various numerical techniques. Examples of practical prob- 
lems are given for which the numerical techniques have been applied, 


NOMENCLATURE 

J 

number of transverse grid points 

A 

matrix, Eq. (54) 

k 

number of time steps 

a k 

cell coefficient, Eq. (53) 

L* 

length of duct, m 

b k 

cell coefficient, Eq. (53) 

m 

spinning mode number 

C k 

cell coefficient, Eq. (53) 

N 

total number of finite element nodes 

c* 

ambient speed of sound, m/s 

N Grid 

numbe f grid points 

C 2 

1 - (Y - l)/2 + 4>f] 


element interpolation function for node n 

d k 

cell coefficient, Eq. (53) 

Nsub 

number of storage locations in sub matrix 


duct diameter 

n 

node number in finite element analysis 

E 

total number of elements 

nt 

transverse mode number 

e 

element 

P* 

pressure, N/m 2 

a k 

cell coefficient; Eq. (53) 

P 

time dependent acoustic pressure, P*/p 0 C 0 

F 

initial condition vector, Eq. (54) 


acoustic pressure, P*(t,x,r), Eq. (12), N/i 

f* 

frequency, Hz 

K 

mean flow pressure, N/m^ 

f 

function of y, see Eq. (39) 

p i 

fluctuating acoustic pressure, N/m^ 

H* 

duct height, m 

p 

spatially dependent acoustic pressure 

I 

number of axial grid points 

Po 

wave envelope pressure, Eq. (62) 

i 

vq: 

R* 

gas constant, m 2 /sec 2 K 



r 

radius dimensionless, r/r * 



* 

*0 

radius "of dUct 



Ar 

radial grid spacing 


T* temperature, K 

T* mean field temperature, K 

T* fluctuation acoustic temperature, K 

t dimensionless time, t*/tp 

tp period lit*, sec 

At time step 

U time dependent axial dimensionless acoustic 

velocity, U(x,y,t), U*/c£ 

U 0 dimensionless mean flow velocity in axial 
direction, U*/C* 

u spatial dimensionless axial acoustic velocity 
u(x,y) 

V time dependent dimensionless transverse or 
radial dimensionless acoustic velocity, 
V(x,y,t) 

-f 

V* vector flow velocity, m/suc 

V* mean vector flow velocity, m/sec 

V* fluctuating vector flow velocity, m/sec 

v spatial dimensionless transverse or radial 

dimensionless acoustic velocity, 

W time dependent dimensionless acoustic velocity 

in 0 direction, W(x,r,0,t) 
w spatial dimensionless acoustic velocity in 
direction, w(x,r,0) 

x axial coordinate, x*/H* or x*/r* 

x m mapped axial coordinates 

Ax axial grid spacing 

A ★ 

y dimensionless transverse coordinate, y /H 

y m mapped transverse coordinate 

Ay transverse grid spacing 

Z impedance, kg/m sec 

Y ratio of specific heats 

C specific acoustic impedance 

ri dimensionless frequency, H*/f*/C* (cartesian) 

or r*f*/C* (cylindrical) 

X dimensionless wave length 

0 angle, radians 


6 r dimensionless specific resistance 

p dimensionless acoustic density fluctuation ,p*/p* 

p* density, kg/m 3 

P* ambient air density, kg/m 3 

p* fluctuation acoustic density, kg/m'* 

flow potential function, m 3 /sec 
mean flow potential function, m 3 /sec 
<f' dimensionless acoustic potential function, 

$*/C*dJ 

spatial acoustic potential function, m 3 /sec 
<f>* fluctuation velocity potential function, 

n 

nr/sec 

X dimensionless reactance 

oj angular frequency, w*d*/C* 

w* angular frequency 

SUBSCRIPTS 

e exit condition 

i axial index (fig. 2) 

j transverse index (fig. 2) 

0 ambient conditfcm or mean flow variables 

s spatial value 

t transient 

1 fluctuation quantity 

SUPERSCRIPTS 

* dimensional quantity 

k time step index 

(1) real part 

(2) imaginary part 

C) vector quantity 

introduction 

With the introduction of strict aircraft noise 
regulations in the late 1960's, the new aircraft 
nacelle designs required acoustvc treatment in the 
inlet and exhaust ducts to reduce engine fan noise. 

To minimize the weight penalty of wall treatment, the 
aerospace industry has been concerned with reducing 


2 


che lengch of a liner for a required sound attenua- 
tion. To perform this optimization, sophisticated 
and flexible suppressor analytical techniques are 
needed to handle sound propagation in ducts with 
axial variations in cross-sectional area, wall liner 
impedance (absorbers), and with gradients in the flow 
Mach number. In an attempt to meet this need, both 
finite difference and finite element techniques were 
developed. In a parallel effort, universities and 
industry have also developed similar numerical pro- 
grams to analyze acoustic propagation in mufflers and 
ventilation systems. 

In this paper, a review is presented covering 
both finite difference and finite element analysis of 
small amplitude (linear) sound propagation in straight 
and variable area ducts with and without flow. Both 
"steady" state and transient theory will be discussed. 
In the "steady" state theory, the pressure and acous- 
tic velocities are assumed to be simple harmonic func- 
tions of time; thus, the governing linearized gas- 
dynamic equations become independent of time. 

The numerical analysis to be considered employs 
a full two or three dimensional finite difference or 
finite element theory. The theories to be discussed 
will have provision for inlet and outlet flow, such 
as in a muffler or jet engine duct. 

Application of finite difference or finite ele- 
ments to acoustic propagation in enclosed rooms 
(refs. 1-8) or to coupled structural-acoustic radia- 
tion systems (refs. 9-12) will not be included. Also, 
theories will not be covered which use the finite 
difference analysis to determine the eigenvalue and 
corresponding modal pressure difference across the 
duct (refs. 13-17). These works have been included 
in the list of reference for completness. 

First, the equations governing sound propaga- 
tion are presented along with the boundary and ini- 
tial conditions associated with a typical turbojet 
engine duct. Next, the various numerical theories 
developed in the literature are presented in Che 
following order: 

(1) Steady State Finite Difference Theory 

(2) Steady State Finite Element Theory 

(3) Special Numerical Transformations 

(4) Transient Numerical Theory 

Items (3) and (4) represent special attempts to nllm- 
inate the problems associated with the highly oscil- 
latory nature of sound propagation in the inlet of a 
turbojet engine, 


and neglecting viscosity and heat transfer, these 
equations are: 


Continuity 

f£+ v*' (p*V*) - 0 

(1) 

Momentum 

p * ^|“y + V* • V^*) - - y*P* 

(2) 

Energy 

& + ■ jh ($ fc *) 

(3) 

State 

P*/p* - R*T* 

(4) 


Where an asterisk denotes dimensional quantities. 


Linearized Gas-Dvnamlc Equations 

Customarily, the dependent flow parameters p*, 

V , and F* are represented as the sum of a mean 
quantity and a small fluctuating quantity, i.e. 
p* - pJ 5 + p*. t* * T q + Tf , $* - V* + Vf and 
P* “ V* + P* where the subscript o represents the 
mean component and 1 represents the fluctuating 
component. The mean flow quantities are assured in- 
dependent of time and thereby satisfy the steady mean 
flow equation. Substituting these assumed expressions 
into Eqs. (1) to (4), eliminating the mean flow quan- 
tities and neglecting nonlinear acoustic quantities 
yields: 

. * 

ipi * *** 

Continuity + V* • (p*v£ + pjv*) - 0 (5) 


Momentum 
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W + v 0 


V*V* + Vi 


V*V$ 


Energy 


State 


+ pjv* • 7*V* - - V*pJ 


C* 2 (* p * + v* 

C ° \at* + V ° 




Pi 


+ Vi 


-A *1 

7 Pn 


+ cf V* . V*p* 
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Pi 

Po 


'* * * -* it * 

"at* + v* VP* + v x * • VP* 


To* 


( 6 ) 


(7) 

( 8 ) 


GOVERNING EQUATIONS 

In this section, the general fluid flow equa- 
tions will first be introduced in vector form. From 
these equations, the general linearized gas-dynamic 
equations, which govern the propagation of sound down 
a duct, will be developed. These equations will be 
specialized to the case of a parallel shear flow in 
a duct with uniform mean temperature and pressure. 
Next, the parallel shear flow equations will be 
written for cylindrical and cartesian coordinates. 
Finally, the equations will also be presented in 
terms of the velocity potential. 


Unsteady Invisld Flow 

The basic equations governing sound flow in air 
ducts are the usual continuity, inviscid momentum 
(Euler), energy, and equation of state. In the 
absence of any external sources of mass or momentum, 


The vector form of Eqs. (5) to (7) were listed by 
Goldstein (ref. 18, p. 5). These equations are fre- 
quently referred to as the linearized gas-dynamic 
equations . 

For the special and important case where the 
mean velocity remains constant and parallel to the 
duct walls (parallel shear flow) and for which the 
mean flow temperature and pressure are constants, the 
energy Eq. (7) reduces to: 

- C fpj (9) 


This derivation is most easily accomplished by first 
expressing the energy equation in terms of entropy 
(ref. 18, p. 8). The linearized continuity and momen- 
tum equations become: 


3P 1 *2 a * 

at* + c o p o 7 
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Mix 

S o \3t* 


+ vS 


V V£ + VJ 


s) - - V*P* (U) 
(Parallel shea; flow) 


*' 

vv 


Cylindrical Coordinates 

Figure 1 shows the typical ducting for a turoo- 
fan engine. Noise from the fan travels up the circu- 
lar inlet against the flow and down through the 
annular exhaust duct as shown. Thus, a cylindrical 
coordinate system (x,r,0) is appropriate for the 
analysis. In general, the acoustic field is three 
dimensional. The sound variation in the axial x 
direction will depend on frequency, duct area, and 
be very large in the presence of wall absorbers. The 
acoustic field may also vary radially from the hub to 
tip at the rotor (see fig. 1), Because of the rota- 
tion of the rotor blades, large circumferential vari- 
ations (see fig. 1) in acoustic pressure may also 
occur depending on blade number and engine rpm. 

A three dimensional solution for sound propaga- 
tion would be expensive to perform. Customarily, 
since the equations are linear, the circumferential 
acoustic pressure and velocity variations are decom- 
posed into spinning modes m: 


Pj(x,r,6,t) - ]£ P^x.r.tje 1 * 6 (12) 

m 

The summation is over those lobe numbers (m) that are 
likely to occur in a particular application. Using 
the assumption that 

P* and v* » e imB (13) 


Equations (10) and (11) for parallel shear flow in 
cylindrical coordinates become after nondimensional- 
izing 
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Equations (14) to (17) are in dimensionless form, 
where the particular nondimensionalizing variables 
are tabulated in the list of symbols. The dimension- 
less frequency n based on the radius of the duct is 
given by 


n 


mil* 

2tt C* 


/ fV > 

lor — jf- for cartesian coordinates to be used later 
\ C 0 t 


(18) 


Equations (14) to (17) will combine to yield a 
third order partial differential equation in pressure, 
as shown by Savkar (ref. 19). In the absence of shear 


("gjjp * O'), Eqs, (14) to (17) combine into the usual 
plug flow wave equation 
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2 nUo (19) 

l dp 3t3x 




In the Bteady state analysis, the noise source 
and dependent acoustic parameters are always assumed 
to be of the form 


p* or V*-e iu * t * « ti 12nt 


(20) 


Consequently, Eqs. (14) to (17) and Eq* (19) reduce 
to 


-g-g-£ w - f -„ 0 jj (2D 

(U^u.-g-i^g-^v (22) 

(12nn)v - - g - U 0 g (23) 

(i2nn)w - - iE p - U 0 (24) 

r 3x 


(2rr,) 2 P - C - U 2 ) p - 14nnU 0 |E (25) 

0 3x 2 3T 2 r 2 t 0 3x 

4- — • *-*“ 

r dr 

The acoustic variables shown in lower case represent 
the time independent values. 


Cartesian Coordinates 

The annular exit duct in the turbojet engine can 
be represented approximately by rectangular geometry. 
Also, many experiments are performed in rectangular 
geometries, For simplicity, spinning modes have been 
neglected. In the absence of spinning modes, the rec- 
tangular forms of the governing equations are 

3P 1 JU 13V HalL . . 

at “ ” IT ax " q 3y - n 3x 
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3V _ _13P_iJo_3V 
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(28) 


Likewise, the steady state equations (using eq. (20)) 
become 
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For the. case of no shear (jTiy’ “ o), Bqs. (26) 
to (28) combine into a single second order wave 
equation 
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The steady state form of Eq. (32) becomes 

(1 - U 2 ) iHS + - iA«r|U 0 E. + (2nn) 2 P - 0 

3x z 3y z 3x 


(32) 


(33) 


l(*) - (C 2 - 4.J)4- rr + (C 2 

+ [- (Y + D*rr»r - 

- (Y - D« r «xxJ *r + 
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- 1) I 4- r * x J$ x 
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In the absence of flow (U 0 equals zero) Eq. (33) re- 
duces to the classic Helmholtz equation 


+ iw(Y - 1) 
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(38) 
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•Mj q. — —P + (2sn) 2 p * 0 (34) 

3x z 3y z 


Potential Flow Equations 

When using numerical techniques in calculating 
sound propagating in ducts, the velocity potential 
formulation of the acoustic wave equation offers many 
advantages over the conventional linearized gas equa- 
tions. For three dimensional flows, the velocity 
potential approach reduces the computer storage and 
running time by an order of magnitude compared to the 
more general linearlized gas equation approach, 
because only one dependent variable $ is required. 
Since the flow into the inlet is usually modelled by 
potential flow (excluding the boundary layer) , the 
acoustic velocity potential is ideally suited for 
acoustic inlet calculations. Sigman, Majjigl and 
Zinn (ref. 20) give the potential gas-dynamic equa- 
tions for a cylindrical duct. The general velocity 
potential equation for inviscid flow (ref. 16, 
eq, (1)) 


4. JL (7*5* . 7 *j*) 
at: *2 'at* 1 9 ; 


+ 1 y*5* • V*(vV . vV) -c* 2 V*V* . o (35) 


was the starting point for their analysis, First, 
they assumed 

$*■♦* + ifij (36) 

where again the symbol t* represents the steady 
mean value of the velocity potential and the 

fluctuation quantity. Next, because the fan gener- 
ates spinning acoustic modes, they assumed 

**(x,r,6,t) - ^(x,r)e- iut e im0 (37) 

Substituting Eqs, (36) and (37) into Eq. (35), yields 
(dimensionless) 


BOUNDARY CONDITIONS 

To obtain a solution to the previously developed 
equations, boundary conditions are required at the 
liner entrance, duct liner walls, and exit plane. 


Entrance Condition 

For a rectangular duct, the boundary condition 
at the liner entrance P(o,y,t) can be of form 

P(o,y,t) « f(y) e L*t* (39) 

where f(y) represents the transverse variation in 
pressure. Similar expressions for $ or acoustic 
velocity can be used. In cylindrical coordinates, 
the radial coordinate r would replace y in Eq. 

(39). In a rectangular duct, spatial pressure vari- 
ations at the liner entrance f(y) have been taken to 
be plane (f * 1) or some harmonic (cos[ntry]) associ- 
ated with the normal-mode analytical solutions (ref. 

21, p. 504), In cylindrical ducts, the spatial vari- 
ations in the radial direction are of the form of 
Bessel functions (ref. 21, p. 511), 

Physically, boundary condition (39) represents 
the sum of a forward and reflected acoustic wave at 
x » o, In a uniquely different approach, Eversman 
et al. (ref. 22) assume an infinite duct upstream of 
the source, and that Eq. (39) represents only the for- 
ward propagating wave. Their analysis separates the 
reflected and transmitted waves. At the present time, 
both approaches to the entrance conditions are used in 
modelling the sound input into the acoustic suppressor, 


Wall Boundary Conditions 

The boundary condition at the surface of a 
locally reacting sound absorbant soft-wall duct can be 
expressed in terms of a specific acoustic impedance 
defined as the ratio of pressure to transverse acous- 
tic velocity: 
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JL. . 1 „£ 

p*C* V v 
o o 


(40) 


For shear flow, the mean velocity at the wall is zero; 
therefore, substituting Eq. (40) into Eq. (28) yields 


For harmonic excitation, e 


ii/t*. 


Eq. (41) reduces to 


la M . i2aa p 

3y c p 


(42) 


Equations (41) and (42) are useful when working with 
the wove equation formulated in terms of only pres- 
sure. Equation (40) is more useful when solving the 
continuity and momentum equations simultaneously. 

It is also convenient to express the specific 
acoustic Impedance in terms of resistance 0 r and 
reactance x as 


t - 9r + iX 


(43) 


Exit Impedance 

In a manner similar to the wall impedance, the 
axial impedance at the duct exit can be defined as 


Ce 


Ua*/H*,y,t) 


(44) 


General values for the exit impedance for single 
modes con be found in references 20 or 23, As with 
the source condition, Eversman, et al. (**£- 22) use 
the technique of accounting for reflected and trans- 
mitted waves. For arbitrary multi-modal wave forms, 
however, a general impedance equation is not avail- 
able for use in the numerical techniques. 


Initial Conditions 

In the transient analysis, for times equal or 
less than zero, the duct is assumed quiescent, that 
is, the acoustic pressure and velocities are taken to 
be zero. For times greater than zero, the applica- 
tion of the noise source (eq, (39)) will drive the 
pressures in the duct. 


STEADY STATE FINITE DIFFERENCE THEORY 

In 1973, Alfredson (ref. 24) and Buumelster and 
Bittner (ref. 25) presented a finite difference theory 
for noise propagation in both hard and soft wall ducts 
in the absence of flow. In the 1973 Aero-Acoustic 
conference, Baumeister and Rice (ref. 26) simulated 
noise propagation in the jet engine duct by Including 
mean flow with the finite difference theory. 

In the first step of these finite difference 
analyses, the continuous acoustic flow field is lumped 
into a series of grid points for which discrete values 
of the acoustic pressure are defined, as shown in 
Fig. 2 for a rectangular duct. 


Grid Point Requirements 

The question naturally arises as to how many 
grid points are required in the axial (x) direction 
and the transverse (y) direction to obtain reasonably 
accurate results. The following one-dimensional ana- 
lysis will be helpful in determining the required 
number of grid points, For the special case of plane 
wave propagation in a hard wall duct, no pressure 
gradients will exist in the transverse direction in 
the duct; consequently, the wave Eq. (33) reduces to 
the one dimensional form 


For a semi-infinite duct (no reflections), the solu- 
tion to Eq. (45) for the complex acoustic pressure is 


P “ 



(46) 


Thus, the pressure oscillates down the duct with a 
wave length 


x.ilia 

* n 


(47) 


The real part of the acoustic pressure in 
Eq. (46) is shown in Fig. 3 for the cije of a dimen- 
'.iortless frequency n “ 5 and a zero Mach number. 
Five separate oscillations appear in the duct. If a 
Mach number U 0 of -0.5 is introduced to simulate 
the inlet flow of a turbofan engine, Eq. (47) indi- 
cates that the wave length would decrease by half, 
and ten full oscillations would appear in the duct. 

In a detailed study of the accuracy of the numerical 
solutions to this problem, Baumeister and Bittner 
(ref. 25) showed that at least twelve grid points per 
wavelength are required to accurately resolve the 
pressure and the transmitted acoustic power. Thus, 
the number of grid points required in the axial di- 
rection is 


I > 13 L*/H* n 
“ 1 + U 0 


(48) 


Besides the axial acoustic oscillations, trans- 
ferse acoustic oscillations (ref, 21, p. 492) con also 
exist in the duct. In a recent paper, Baumeister 
(ref. 27) showed that the number of grid points J in 
the transverse direction required to accurately re- 
solve the highest order propagating mode is 

J > 12n (49) 

Thus, the total number of grid points N Gri(J required 
in the duct is the product of Eqs. (48) and (49): 


N Grid 


144 n 2 L*/H* 
1 + U 0 


(50) 


Although a large number of grid points is re- 
quired particularly at high frequence s, a more diffi- 
cult problem results from the natur . of the governing 
steady state acoustic equations themselves, as will be 
considered next. 


Difference Equation 

The continuous governing differential equations 
can be changed to a system of algebraic equations by 
means of the finite difference approximations. This 
set of algebraic equations can then be solved simul- 
taneously to determine the pressure at each grid point. 
In constructing the difference equations, only sound 
propagation in the absence of flow will be developed 
here in order to keep the mathematics simple, 

Without flow, the classic Helmholtz equation 
(eq. (34)), governs pressure propagation in the duct. 
Using the usual 5 point difference approximations, 

Eq. (34) becomes 

A + Pi+l.ij + ^Pl.i-1 ~ 2 Pl.1 + Pl.i+lj 


( ! 


Pj-l.i " 2 Pl 


Ax 


(45) 


+ (2trn) p it j 
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(1 - U 2 ) ^ - i4rnU 0 + (2Ttn) 2 p = 0 
dX 


0 


(51) 


or for Ax equal to Ay 

p i-l,j +p i, j-1 - B - C2*n4x)Z]p t(j + P 1(j+l + P i+ltJ - o 

(52) 

Equation (52) applies only to cell 1 in Pig. 2 which 
is not adjacent to the boundaries. For the other 
cells in Fig. 2, a complex impedance condition is 
specified along a boundary. The integration method 
is most convenient to use for generating the differ- 
ence equations at the boundaries, as fully documented 
in Ref, 28, Appendix D, 

The finite-difference approximations for the 
various calls shown in Fig, 2 con be expressed in 
terms of the coefficients 

Wl,J + Vi.J-l + c K p i.j + Vi,j+l + Vl+l,j"° 

(53) 

The subscript k denotes the cell number. The col- 
lection of the various difference equations at each 
grid point forms a set of simultaneous equations that 
can be expressed as 

(A}[P] » [F] (54) 

where (A) is the known coefficient matrix, fjH is the 
unknown pressure vector, and U;J is the knoVn column 
vector containing the various Boundary conditions. 

The matrix is complex because of the complex nature 
of the source and impedance boundary conditions. 


Matr ix Solution 

The magnitude of the frequency term (2«nAx) 2 in 
Eq, (52), which subtracts from the main diagonal ele- 
ment ck of coefficient matrix represented by 
Eq, (54), is such that the matrix is not diagonally 
dominant. As a result, Iteration solutions cannot be 
used. Equation (54) is usually solved by the Gauss 
elimination technique, In general, the elimination 
technique requires the storage of all the matrix 
elements which would be represented by twice (complex 
nature) the square of Ng rll) given by Eq. (50). For 
high frequencies (n) , this represents a large burden 
on tne storage capacity of a computer. Sparse-matrix 
techniques (ref. 29) have been employed to reduce 
computer storage and run times as much as possible, 
Quinn (ref, 29) partitioned the matrix (A) into a 
tridiagonal form, which reduced the elements in the 
submatrix to 

N sub - J 2 - 144 n 2 (55) 

Quinn's program stores the elements of the sub- 
matrices in separate locations on tape or disk. 

These matrices are then read in one at a time to ob- 
tain a solution of the matrix. 

Baumeister (ref, 27) and Quinn (ref. 30) have 
modified similar matrixes to allow iteration tech- 
niques; unfortunately, the convergence is too slow to 
be of any practical value, Other approaches, such as 
in Ref. 31, might still offer iterative possibilities. 


Example Solutions 

Figure 4 shows some sample numerical calcula- 
tions for a plane wave propagating in a uniform hard 


wall duct, As seen In Fig, 4, the numerical theory 
and exact analysis are in good agreement. 


AmittiM t 1 m. 

So far in this paper, the difference theory con- 
sidered only straight uniform area ducts. Generally, 
finite differences when applied to irregular ducts 
requires considerable bookkeeping (ref, 32, p, 363). 
Quinn (refs, 33 and 34), however, applied the method 
of conformal mapping in conjunction with finite dif- 
ferences to analyze sound propagation in a hyperbolic 
horn shown in Fig, 5. In the absence of flow, the 
following transformations 

x . • Tf-S cos y m (56) 


a*” 1 + e"*® 

y „ — - sin y m (57) 


change the governing Helmholtz equation into 

■M 1 + + (2itn) 2 F(x. ,y )p » 0 (58) 

K K 

where 

e 2 *® - e” 2 *® + 2(coe 2 y m - Sin 2 y m ) 

p (*i#»yir,) “ 4 — — — (59) 


Equation (58) was programmed in finite difference 
form hy Quinn for the straight duct shown in Fig. 5, 
Quinn also conformally mapped a conical horn and 
found the numerical results to be in excellent agree- 
ment with known analytical solutions. 

In general, for ducts with variable area, the 
finite element theory ia the most convenient to use. 
However, ab pointed out by Quinn (ref, 30), the advan- 
tage of finite elements versus finite differences Is 
not as great when a mapping function is used to com- 
pute the mean flow field. For example, conformal 
mapping has been used In the calculation of subsonic 
and transonic inlet flow fields in the DC-8 engine 
(ref. 35). In such a case, where the conformal trans- 
formation is available, the use of finite difference 
theory would clearly have an advantage over finite 
element theory. 


STEADY STATE FINITE ELEMENT THEORY 

The use of finite element techniques in duct 
acoustics began in 1975. Young and Crocker (ref, 36) 
developed a two-dimensional finite element solution 
of acoustic propagation in mufflers without flow using 
a variational technique. In a preliminary investiga- 
tion, h '.pur and Mungur (ref. 37) concluded that a 
finite element Galerkin method could be successfully 
employed to handle shear flow in ducts with area var- 
iations. Since these papers, numerous finite element 
solutions of the duct propagation problem have appear- 
ed in the literature and will be discussed herein. 
First, however, a brief review of the finite element 
theory is presented to establish categories to clas- 
sify the various references in the acoustic literature. 
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Basic Theory 

The solution of the duct acoustic problem by 
the finite element theory follows the usual step-by- 
step process (ref. 38, p. 7) developed for any con- 
tinuum problem. These step* will now be succinctly 
listed with particular emphasis on how they relate 
to the duct acoustic problem. 

1. Discretize the Continuum: The first step 
is to divide the interior of the duct into elements. 

A variety of elements have been uB«d in the duct 
acoustic problem as shown in the first column of 
Table I, Table I will be used os a focal point about 
which the discussion of the acoustic literature will 
center, 

2. Select Interpolation Function: The next 
step is to assign nodes to the element and then 
chooBe the type of interpolation function to repre- 
sent the variation inside the elements of the field 
variable, such as pressure and the acoustic veloc- 
ities. In Table I, column two shows the node pattern 
while column three describes the interpolation func- 
tion. The required number of elements can be estim- 
ated from the results of the finite difference theory. 
For sound propagation without large reflections (ex- 
clude low frequency as in mufflers), Eq. (SO) con be 
used to estimate the number of finite elements neces- 
sary to obtain an accurate solution. In determining 
the required number of finite elements, however, the 
constant 144 will be smaller because the interpola- 
tion function in finite element theory is generally 
of higher order than the first order linear differ- 
ence approximation used to establish Eq. (50). The 
precise decrease will depend on the type of element 
used in the analysis. In general, the more nodes per 
element the fewer will be the number of elements 
needed to resolve the acoustic field. On the other 
hand, since the size of the solution (global) matrix 
is proportional to the number of nodes, considerably 
more computational time and computer core memory are 
required for the higher order elements. Column four 
lists the number of dependent variables considered 

in the particular study. For a fixed number of ele- 
ments, the computer storage will be proportional to 
the square of the number of dependent variables while 
the solution times can increase by the cube (ref. 39, 
p. 261) of the dependent variables. Clearly, those 
equations should be employed which can successfully 
model the physics with the fewest dependent variables. 

Generally, Lagrange polynomials are used for 
interpolation when only one field variable (such as 
pressure p) is calculated at each node, commonly 
called the C° continuity problem. When using Hermite 
polynomials, the values of 3p/3x and 3p/3y can 
also be calculated at each node, commonly called the 
Cl problem. In general, cl continuity is more diffi- 
cult to construct (ref. 38, p. 172), so with the ex- 
ception of items 1 and 5 in Table I, Lagrange poly- 
nomials are used. Also, if Cl continuity increases 
the unknowns by four (item 5, Table I), then the 
resulting solution matrix will increase by about 
sixteen. 

3. Finite Element Method: Next, the unknown 

nodal values of the pressure (or any other dependent 
variable) must now be constrained such that the nodal 
pressures provide an accurate approximation for the 
true pressure distribution. In conjunction with the 
governing differential equations and boundary condi- 
tions, a finite element method is used to develop the 


algebraic element equations lot the unknown pressures 
associated with each element. The element equations 
are then combined into a general matrix (global) which 
is solved for the pressures at the nodes. 

In acoustics, except for the no flow case, the 
Calerkin finite element method is most often used to 
establish the matrix equations for the unknown field 
variable, For example, consider the case using the 
velocity potential to describe the acoustic field, 

Eq, (38). Application of the Calerkin Finite element 
method yields the following relationship 

E 

jj N®L(*)dA< a > - 0 n - 1.2....N 


where the integration is performed over the area of 
each element. The known interpolation function N£ 
is for the n node in the e th element and N is 
the total number of nodes in the problem under consid- 
eration, It should be noted that N® is zero for all 
elements not having the nodal point as a vertex. 
Equation (60) provides N equations for the N 
unknown nodal values. 

4. Matrix Equations' To find the solution to 
the complete network of elements, all the element 
equations must now be assembled into a global matrix. 
With the exception of the least squares approach, the 
solutions are all accomplished by soma sort of elim- 
ination procedure. Effective solutions take into 
account the banded nature of the matrix. The band 
width increases with a larger number of element nodes 
and number of dependent variables. As with finite 
dixterence, the storage requirement of the finite 
element solution can severely restrict the range of 
problem that can be solved as pointed out by Quinn 
(ref, 30), For the larger problems, an out-of-core 
banded solver is generally required with a moderate 
amount of in-core storage but much more input/output 
time, 


Finite Element Theories 

The following discussions of the published fin- 
ite element theories will be in order of the complex- 
ity of the flow field and interpolation function and 
not chronological. 

No Flow , Young and Crocker (ref. 36, item 1, 
Table I) developed a two-dimensional finite element 
model for analyzing sound propagation in expansion 
chambers with hard walls. The basic elements were 
rectangles using one-dimensional Hermitian polynom- 
ials of fourth .order. Good agreement was found 
between the finite-element approach and the acoustic 
filter theory at low frequencies where plane wave 
analysis applies. At the higher frequencies, where 
the diameter to wave length ratio is greater than 0.8, 
significant differences between the plane wave and the 
two-dimensional finite element analysis were found. 

Kagama and Qmote (ref. 40, item 2, Table I) ex- 
panded the analytical work of Young and Crocker to 
lined axisymmetric soft wall mufflers with more com- 
plicated inputs. A linear three node triangular 
element was used in Ref. 40 for the hard wall muffler 
while a six node element triangle was used in Ref. 41, 
item 3, Table I, in the soft wall analysis. A wealth 
of example problems were considered in these papers. 
Rectangular and conical mufflers were analyzed with 
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and without sound absorbing materials. Significantly, 
experimental data were obtained for all the examples 
analyzed. Excellent agreement between the finite 
element theory and experiments was obtained. There- 
fore. Che acoustic response of a muffler of any shape 
can be accurately and easily determined by a two- 
dimensional finite element analysis. In the same 
ti«ce frame, Craggs (refs. 42 and 43) presented a 
formulation and approach similar to that of Kagawe 
and Omute. Craggs used a more refined axlsymraetrlc 
hexahedral isoparametric element as shown in Table I, 
item 4, Isoparametric elements allows the sides of 
the hexahedral element to bend and more closely match 
the shape of a curved boundary. In Refs. 44 and 45. 
Craggs extended his treatment of sound absorbing 
materials to include extended reaction effects. His 
two-dimensional model showed that propagation parallel 
to the incident boundary can be important 

Tag, Lumsdnine, and Akin (refs. 46 and 47) have 
applied linear rectangular isoparametric elements to 
two-dimensional (cartesian) rectangular ducts with 
soft walls without flow. Examples are shown 
(ref, 46) for sound propagation around circular and 
rectangular bends as might commonly be seen in indus- 
trial ducts, Some simple plane wave propagation 
examples with plug flow in hard walls (ref. 46) are 
also presented. Similar to Young and Crocker, Lester 
and Parott (ref. 48) have recently presented a 
Hermitian analysis for sound propagation In straight 
ducts with sound absorbing walls. Finite elements, 
Wiener Hopf, and wave propagation math models were 
found to glva consiutant results. 

Irrotatlonal Flow . Sigmon* Majjigi and Zinn 
(ref. 20, item 7) presented the first finite element 
analysis which treated irrotatlonal flows in ducts, 
They used a simple linear three node triangular 
element in conjunction with Galerkin analysis to 
model the sound propagation. The potential function 
wave equation is appropriate to inlet flow in a tur- 
bojet engine where the flow is nearly irrotatlonal 
(excluding the boundary layer). The potential flow 
formulation gives an order of magnitude reduction in 
required core storage compared to the more general 
flow solutions to be presented next. As an illustra- 
tive example in their paper, they analyzed the hard 
wall version of the NASA QCSEE (quiet, clean short- 
haul experimental engine) inlet. This inlet was de- 
signed to suppress inlet-emitted engine machinery 
noise using a high throat Mach number and thus has 
somewhat more curvature than a conventional inlet. 

For an initial plane wave, their finite element 
analysis showed that significant distortion of the 
initial plane wave occurred due to the entrance curv- 
ature effects. 

Majjigi, Sigman, and Zinn (refs, 49 and 50, 
item 8) have expanded their original analysis to 
include sound absorbing material by using a quadratic 
interpolation function. Using this program, Baumeis- 
ter and Majjigi (ref. 51, item 8) estimated the sen- 
sitivity of duct attenuation to inlet curvature, 
centerbody and irrotational flow gradients for a high 
Mach number turbojet inlet. They showed that area 
variations can significiantly change the calculated 
attenuation as compared to the circular cylinder 
model presently being used for the design of such 
inlets. 


The assumption of lrrotationality eliminates 
tho practical case of parallel shear boundary layer 
flow, Baumeistar and Majjigi (ref. 51) showed that 
even in an approximate manner the velocity potential 
cannot be used to estimate the effects of sheared 
flow, However, a boundary layer correction based on 
the Goldstein-Rice analysis (ref. 52) was suggested 
for predicting the attenuation with shear in a soft 
wail duct. 

Tag and Lumsdnine (ref, 53, item 5) also per- 
formed a similar irrotatlonal analysis for a hard wall 
inlet using an isoparametric quadrilateral element. 
They showed that the finite element analysis can con- 
veniently handle the problem of cutoff in a variable 
area duct. For a sample inlet, they also performed 
calculation* to show the effect? of duct geometry on 
far field attenuation and of flow gradients on the 
effective sound wavelength. 

General Flow Field , Items 10, 11, and 12 in 
Table I are concerned with formulating a general fin- 
ite element solution of the complete acoustic flow 
field, Eqs. (5-8). In most cases, the general theory 
hos been applied to some relatively simple cases to 
check the theory. 

Quinn (ref. 30, item 10) examined the use of 
various interpolation functions and finite element 
methods on many sample problems. He compared the 
accuracy of linear interpolation functions on tri- 
angles, bilinear interpolation function on rectangles 
and biquadratic interpolation functions on rectangles, 
He found that the biquadratics permiu good approxima- 
tion of curved boundaries and better convergence than 
the bilinear interpolation functions. 

Abrahamson (refs, 54-56, item li) has developed 
a finite element model using a linear rectangular ele- 
ment along with the Galerkin method . He has applied 
the method to conical horns, liners with axial imped- 
ance variations and the optimization of liners for 
maximum acoustic attenuation. He has devoted consid- 
erable effort (ref* 56) to reducing the computational 
times for large complex problems. 

Ascley, Evarsman, and Thanh (refs. 57-60, 
item 12) have developed a finite element method based 
on weighted residuals using an eight node isoparamet- 
ric element. As in Ref. 30, they also found the 
Galerkin method to be superior to the least squares 
approach. For uniform rectangular and circular ducts, 
they showed that the finite element method produces 
results for eigenvalues and transmission coefficients 
nearly identical to those from analytical approaches. 
Eversman, Astely and Thanh (ref. 57) also used the 
finite element method to explore the choking that 
might occur in a subsonic flow. They found that a 
large transmission loss is not an automatic conse- 
quence of propagation against a high subsonic mean 
flow. 

NASA Lewis Research Center has begun an experi- 
mental program to check the results of the finite 
element models against experiment. A 0.5 area con- 
traction in the form of a quartic curve has been 
constructed in which experimental results are cur- 
rently being obtained. Figure 6 shows the comparison 
of the Astley-Eversman finite element program with 
the experimental data. As seen in Fig. 6, experiment 
and theory are in good agreement, The finite differ- 
ence theory (dashed line) by Prof, James White will 
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be d incunKed later in the section on transient solu- Ay > ($ 3 ) 

tlons to the wave equation. n 


SPECIAL NUMERICAL TRANSFORMATIONS 

In the inlet to i turbojet engine, the dimen- 
sionless frequencies n ore of order 30 to 50 for 
the higher harmonics of the blade passing frequency. 
The storage requirement and associated running times 
for these high frequencies are not practical at the 
present timi. Some special techniques hove been 
developed in order to reduce the strrage require- 
ments. 


Wave, Envelope T echnique 

To remove some (if not all) of the axially 
oscillatory part of the wave preseure profile and 
thereby greatly reduce the required number of grid 
points, Baumelster (rcfu. 61 and 62) assumes 


p(x,y) « P 0 (x,y)e" i2 * xA (61) 

where p represents the pressure of the oscillating 
curve in Fig. 7 and p 0 represents the pressure of 
the wave envelope shown by the dashed line. The 
length > represents the effective axial wave len;. h 
of the oscillation, In general, X is chosen to fcc 
the value associated with a particular mode propagat- 
ing through a hard wall duct. For the cose of a 
plane wave propagating in a soft wall duct without 
flow, x is approximated by 1/n, Substituting 
Eq, (61) into the Helmholtz Eq. (34) yields the wave 
envelope equation 


- |2 Po 
3x 2 3y 2 


a Pn 

iAitri » 0 

Sx 


(62) 


The finite difference form of this equation was then 
programmed (ref. 60). Excellent agreement between 
the analytical and wave envelope finite difference 
calculations are shown in Fig. 8, For the n “ 5, 
L*/Ir ■ 6 case in Fig. 8, the conventional finite 
difference theory required 3600 grid points, whereas 
the wave envelope difference theory required only 100 
grid points. Thus, a savings of 3500 grid points 
over the conventional difference theory was obtained. 
In Ref. 63, this technique was used to optimize multi- 
element liners of long lengths at high frequencies. 
This technique can also be used with finite elements. 

At the present time, this technique has been 
applied only to the simple cases of no flow and plug 
floe. For a variable area duct with flow, it will be 
necessary to have X a function of x. Calculated 
trial values of X could be obtained from a one- 
dimensional solution of the governing equations. 


Marching Techniques 

Although initial value numerical solutions of 
elliptic equations are generally unstable (ref. 64, 
p, 63), Baumelster (refs. 65 and 66), using the 
Von Neumann stability analysis, showed that the two- 
dimensional Helmholtz wave equation can be solved 
approximately using explicit marching techniques. 
However, the grid spacing in the transverse y 
direction must be sufficiently large so that non- 
propagating acoustic modes cannot be resolved, that is, 


Compared to standard finite-difference or finite- 
element boundary val < approaches, the numerical 
marching technique is ovdare of magnitude shorter in 
computational time and required computer storage, 
This technique is limited, however, to high frequen- 
cies and to cases where reductions are small, 

The numerical optimum attentuaeions were cal- 
culated using the marching technique and compared to 
the corresponding analytical results similar to 
Fig. 8, For dimension!-’’* frequencies n ■ 1 or 
n ■ 2, the numerical results were in poor agreement, 
while for n * 5 or greater the method worked quite 
well. Consequently, marching could be used in con- 
junction with conventional analysis to hsndie the 
special cases of high frequency and long length. 


.p§.r..g|gli.fayp,Uhs 

To obtain the far field acoustic radiation pat- 
terns and to eliminate tbs necessity of specifying 
the impedance of the duct at its exit, finite elements 
con be extended from the internal portion of a duct 
into the far field, as ahown by Kagawo et al. (ref, 

40, fig. 16). In Ref. 40, an acoustic horn with an 
opening into a semi-infinite space was analyzed by 
introducing a hypothetical spherical boundary with a 
characteristics x 0 c 0 exit impedance termination. 
However, when taking account of acoustic radiation in 
this way, o spherical boundary of large area must be 
considered, which correspondingly results in solving 
a large -matrix equation. To eliminate this large 
matrix requirement, Kagawa et al. (ref. 67) recently 
developed a combination of the finite element (in 
duct) and analytical methods (Green-Tbeorm - far 
field) to analyze sound propagation from an ocouatic 
horn. Experiment and theory were found to be in good 
agreement. 

Bnumaistor and Majjigi (ref. 51) presented a 
matrix partitioning approach that separated the duct 
and far field into two or more regions. The regions 
are coupled by on assumed value of impedance. For the 
no flow case in the for field, the dimensionless 
impedance enn he chosen to be unity at a spherical 
surface, In Ref. 51, large matrix storage reductions 
(factor of 50) were obtained. For cases where complex 
flows occur around the exit, the first partition 
should be in the far field where an analytical value 
of impedance can be reasonably specified, 


TRANSIENT NUMERICAL SOLUTIONS 

In order to eliminate the storage requirements 
associated with the steady state numerical analysis 
of sound propagation in a turbojet engine inlet, 
Baumelster devr loped time dependent numerical solu- 
tions for noise propagation in a two-dimensional duct 
without flow (ref, 27), with parallel shear flow 
(ref. 68), and with axisymmetric plug flow (ref. 69). 
Advantageously, matrix storage requirements are com- 
pletely eliminated In the time dependent analysis. 

The solution vectors for pressure and velocity need 
only be stored. 


Difference Equations 

The analysis begins with a noise source radiating 
into an initially quiescent duct. Next, an explicit 
iteration method obtains the transient as well as the 
"steady" state solution. 
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To illustrate how the transient equations are 
programmed consider the no flow cartesian wave equa- 
tion (eq. (29) with Up ■ 0) in one space dimension! 

,» i2t . 2*£ m 


Away from the boundaries, the second derivatives in 
the wave equation can bo represented by th/i usual 
central differences in time (subscript kl and space 
(subscripts i,j) 



where At and Ax are the time and space mash spac- 
ing, respectively. Solving Eq, (65) for PR+1 
yields 1 


pf 1 - 2 pJ - p}' 1 + (^) 2 (j’ 1+1 - 2p i + P J-l) («) 

The procedure is explicit since all the past values 
of pR are known as the new values of P k+ * are 
computed, The derivation with two space dimensions 
and boundary conditions is fully documented in Refs. 
27, 66, and 69, Since Eq, (66) is a simple algebraic 
equation, storage is only required for the solution 
vectors P?, Pk-lj that is, no matrix storage or man- 
ipulation Is required , This was extremely useful In 
Ref. 68, where a simultaneous solution of the contin- 
uity and momentum equations was employed to obtain 
the acoustic values of P, U, and V, 


graphical Output 

Fortunately, the pressure calculated at each 
Instant of time nead not be saved or displayed. 
Although multiple values of pressure are calculated, 
only the latest values (k,k-l) need be stored. After 
steady Btate is reached (checked numarially) , the 
time-dependent results will be compared to the 
"steady" state results simply by dividing by e^at, 
that is 


P(x.y) - EteoZall (67) 

qXaITC 

In this case where the source is a simple harmonic 
function of time ei2ttt, p represents the Fourier 
transform of P(x,y,t) (ref , 18, p, 11). In this 
manner, the numerical results (refs. 27 and 69) are 
nearly identical to those shown in Fig, 4, The book- 
keeping and graphical output is held to a minimum by 
the use of Eq. (67). The time-dependent analysis is 
considerably faster than the steady analysis, as seen 
in Fig. 9, In general, however, the computational 
time associated with tho steady and transient analysis 
will be compatible. The major advantage of the 
transient analysis remains the elimination of matrix 
storage requirements. Because manipulation of 
matrices is omitted, the time dependent approach is 
much easier to program and debug. 


Time-dependent numerical techniques have been 
applied to one-dimensional sound propagation (ref. 70, 
p, 256), two-dimensional vibration problems (ref, 32, 
p. 452), end the more general problem of compressible 
fluid flow (ref, 71), Numerical stability is a major 
difficulty associated with all these time dependent 
solutions. The new problem associated with two- 
dimensional soft wall acoustic problems is the influ- 
ence of the boundary conditions on stability. The 
manner in which the boundary conditions are posed con 
cause the solution to be absolutely unstable 
(ref, 69), At the present time, the time marching 
technique has been found to be stable for no ilow 
(ref* 27), plug flow (ref. 66), and axis symmetric 
spinning wave acoustics in a cylindrical duct 
(ref, 69). Numerical experimentation will be re- 
quired to determine if the equations will remain 
stable and model the acoustic propagation when axial 
and transverse variations in the mean flow are con- 
sidered. 

At the present time, the transient method ap- 
pears to have one major drawback, The transient 
method does not converge for cut-off modes (ref. 69). 
This has implications as to its use in a variable 
area where modes may become cut-off in the small area 
portion of the duct. This will ultimately have to be 
resolved since cut-off modes will be encountered in 
real problems. In the mean time, the steady state 
numerical techniques con be conveniently employed 
near cut-off. Because of the relatively long axial 
wavelength of a near cut-off mode, the required 
number of grid points (or elements) is relatively 
small. 


Variable Area 

The finite element method is inappropriate for 
the two-dimensional transient problems because the 
solutions are implicit requiring storage of the usual 
large matrix. Therefore, in applying tho transient 
method to problems with oren variations, the method 
of mapping in conjunction with the finite difference 
method could be used. 

In particular, the mapping procedure developed 
by Thompson et al. (ref, 72) provides for the auto- 
matic generation of a general coordinate system with 
coordinate lines coincident with all the boundaries 
of an arbitrary shaped duct. Professor James White 
of the University of Tennessee (Knoxville) under NASA 
grant NAG 3-18, is extending the procedure of Ref, 72 
to the problem of aound propagation in variable area 
duets. In a form similar to Eq, (59), Professor White 
has developed the appropriate transformed acoustic 
equations which are solved In the mapped rectangular 
domain. The results are then transformed back Into 
the physical plane by on inverse mapping procedure. 

Figure 6 shows a comparison of Professor 
White's transient finite analysis with the experimen- 
tal data for the quartic duct discussed earlier in 
conjunction with the Astley-Evarsman finite element 
analysis. As seen in Fig. 6, experiment and theory 
are in good agreement, 


CONCLUDING REMARKS 

The finite difference and finite element theories 
have proven to be useful calculational tools in 
handling acoustic propagation in ducts with discontinu- 
ities in lining impedance, with area variations, and in 
ducts where large variations in the mean flow field 
occur. 
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At the present time, the steady and transient 
numerical theories are capable of handling the design 
of mufflers, expansion chambers, and exhaust ducts of 
turbofan engines where the effective frequency n is 
quite small, In contrast, the steady numerical 
theories are not ns yet capable of analysing n prac- 
tical turbofan inle^ because of the high fraquenciea 
and Mach numbers in the inlet, Special steady state 
transformations and the transient technique have 
potential for handling this problem. 

At the present time, NASA, universities and 
industry are extending the steady and transient num- 
erical techniques to situations Involving high fre- 
quencies, high Mach numbers and geometries with 
large area variations, The transient numerical tech- 
nique is being adapted to handle nonlinear affects 
including acoustic shucks. An additional goal is to 
develop mure realistic inlet and exit conditions to 
model a tvplcal turbofan aircraft nacelle, 
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duct for n * 5, L*/H '' * 1. 











Figure 6. - Comparison of Astley-Eversman finite element 
theory and White finite difference mapping theory to 
experiment (n * 0. 172). 
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